
#ifndef POLYMATRIX_HPP
#define POLYMATRIX_HPP

//
// A Matrix of polynomials
//

#include "Poly.hpp"
#include "Matrix.hpp"

using namespace std;

template <class T> class PolyMatrix : public Matrix < Poly < T > >
{

private:

  // disable these methods
  virtual void pop_back () { }
  virtual void push_back () { }
  // row, col & order
  int m_row, m_col, m_order;
  bool m_IsDecomp;
  ios::fmtflags m_ff;

  //protected:

public:

  // access dimensions
  virtual int row() const { return m_row; }
  virtual int col() const { return m_col; }
  virtual int order() const { return m_order; }
  virtual void IOS_Format ( ios::fmtflags p_ff )
  {
    m_ff = p_ff;
  }

  // constructors
  PolyMatrix () : Matrix < Poly < T > > (), 
		  m_row(0), m_col(0), m_order(0),
		  m_ff(ios::fixed), m_IsDecomp(false) { }

  PolyMatrix (int p_row, int p_col, int p_order ) : 
    Matrix< Poly<T> > (p_row), m_row(p_row), m_col(p_col), m_order(p_order),
    m_ff(ios::fixed), m_IsDecomp(false)
  {
    for (int row=0; row < p_row; ++row )
      (*this)[row].resize(p_col);
  }

  //
  // Initialize with a Vector of Matrices
  //
  PolyMatrix ( const Vector< Matrix<T> >& p_M_Coeff ) : 
    Matrix< Poly<T>(p_M_Coeff.size()) > (p_M_Coeff[0].row(),p_M_Coeff[0].col()),
    m_ff(ios::fixed), m_IsDecomp(false)
  {
    m_order = p_M_Coeff.size();
    if ( m_order < 1 )
      throw (range_error("PolyMatrix::PolyMatrix(Vector<Matrix>) : zero order."));
    m_row = p_M_Coeff[0].row();
    m_col = p_M_Coeff[0].col();
    //Matrix< Poly<T> >::resize(m_row,m_col);
    for ( int k=0; k < m_order; ++k )
      {
	const Matrix<T>& Mtmp = p_M_Coeff[k];
	if ( Mtmp.row() != m_row || Mtmp.col() != m_col )
	  throw (range_error("PolyMatrix::PolyMatrix(Vector<Matrix>): Matrix dimension error."));

	for ( int j=0; j < m_row; ++j )
	  for ( int i=0; i < m_col; ++i )
	    (this*)[row][col].Coeff(Mtmp[j][i],k);
      }
  }
  //
  // evaluate for a given x
  //
  Matrix<T> operator()(T p_x)
  {
    Matix<T> Mret(m_row,m_col);

    for ( int j=0; j < m_row; ++j )
      for ( int i=0; i < m_col; ++i )
	Mret[j][i] = (this*)[row][col].eval(p_x);
    
    return Mret;
  }
  //
  // decomposition - sort of
  //
  void Decomp()
  {
    if ( m_IsDecomp )
      return;

    for ( int j=1; j < m_row; ++j )
      {
	Poly<T> DInv = (*this)[j-1][j-1].Inverse();
	
	for ( int i=0; i < m_col; ++i )
	  }

    m_IsDecomp = true;
  }

  //
  // determinant - characteristic polynomial
  //
  Poly<T> Det()
  {
    if ( ! m_IsDecomp )
      Decomp();

    Poly PRet(Vector<T>(1,1));
    for ( int j=0; j < m_row; ++j )
      PRet *= (*this)[j][j];
    return PRet;
  }

  // test for symmetry
  bool virtual isSymmetric() const
  {
    if ( m_row != m_col )
      return false;

    for (int row=0; row < m_row; ++row )
      {
	for ( int col=0; col < row; ++ col )
	  {
	    if ( (*this)[row][col] != (*this)[col][row] )
	      return false;
	  }
      }
    return true;
  }

  // empty matrix
  PolyMatrix& zero()
  {
    for ( int i = 0; i < m_row; ++i )
      {
	for ( int j = 0; j < m_col; ++j )
	  (*this)[i][j] = 0;
      }
    return *this;
  }
  // identity matrix
  PolyMatrix& identity()
  {
    this->zero();
    for ( int i = 0; i < m_row; ++i )
      {
	(*this)[i][i] = 1;
      }
    return *this;
  }

  // diagonal matrix
  PolyMatrix& diagonal ( const Vector<T>& p_vec )
  {
    if ( p_vec.size() != m_row )
      throw (range_error("PolyMatrix::diagonal(Vector) : dimension error"));
    this->zero();
    for ( int i = 0; i < m_row; ++i )
      {
	(*this)[i][i] = p_vec[i];
      }
    return *this;
  }

  // transpose
  PolyMatrix<T> Tr() const
  {
    // result
    PolyMatrix<T> Mtmp ( m_col , m_row );

    for (unsigned row = 0; row < m_row; ++row )
      {
	for (unsigned col = 0; col < m_col; ++col )
	  {
	    Mtmp[col][row] = (*this)[row][col];
	  }
      }
    return Mtmp;
  }

  void write (ostream *os, int p_prec=5, int p_width=10 ) const
  {
    // set output format
    ios::fmtflags oldopt = os->setf( m_ff );
    os->precision(p_prec);

    for (int i=0; i < this->size(); ++i)
      (*this)[i].write(os,p_prec,p_width);
    *os << endl;

    // reset output format
    os->flags( oldopt );
  }

  void read (istream *is)
  {
    // Vector read
    for (int i=0; i < this->size(); ++i)
      (*this)[i].read(is);
  }
};

template <class T> ostream &operator<< ( ostream &os,PolyMatrix<T> &p_matrix )
{
  p_matrix.write(&os);
  return os;
}
template <class T> ostream &operator<< ( ostream &os,const PolyMatrix<T> &p_matrix )
{
  p_matrix.write(&os);
  return os;
}
template <class T> istream &operator>> ( istream &is,PolyMatrix<T> &p_matrix )
{
  p_matrix.read(&is);
  return is;
}

#endif // POLYMATRIX_HPP
